Propagation of coherent waves in elastically scattering media 
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A general method for calculating statistical properties of speckle patterns of coherent waves prop- 
agating in disordered media is developed. It allows one to calculate speckle pattern correlations 
in space, as well as their sensitivity to external parameters. This method, which is similar to the 
Boltzmann-Langevin approach for the calculation of classical fluctuations, applies for a wide range 
of systems: From cases where the ray propagation is diffusive to the regime where the rays experi- 
ence only small angle scattering. The latter case comprises the regime of directed waves where rays 
propagate ballistically in space while their directions diffuse. We demonstrate the applicability of 
the method by calculating the correlation function of the wave intensity and its sensitivity to the 
' wave frequency and the angle of incidence of the incoming wave. 

, PACS numbers: 72.15.Rn, 73.20.Fz, 73.23.-b 

£ , I. INTRODUCTION 

Jo 

Characterization of statistical properties of coherent waves propagating through an elastically scattering disordered 
medium is relevant for a variety of physical situations, ranging from propagation of electromagnetic waves through 
interstellar space or the atmosphere, seismology, and medical imaging by ultrasound or light, to electron transport in 
disordered conductors. When coherent waves propagate through such media their intensity exhibits random, sample- 
i , specific, fluctuations known as speckles. These fluctuations result from the interference of rays traveling along different 
*^ paths. In this article we study the statistics of speckles. 

The problem can be characterized by several length scales: The propagation distance of the ray through the medium, 
Z , the elastic mean free path, £, which is the typical distance the ray travels between two scattering events, and the 
transport mean free path, £ tr which characterizes the typical distance for backscattering. In the limit of very thin 
sample, Z < £, rays move almost ballistically through the sample, since scattering porbability is small. This regime 
has been extensively studied^. In the opposite limit of a very wide sample, Z 3> lt r , the rays propagate diffusively 
in the system. This regime has been considered in Refs. HHGI At spatial scales exceeding the transport mean free 
path the statistical properties of speckles in the diffusive regime (excluding features associated with rare events) are 
characterized by the diffusion coefficient and are independent of the details of the disorder. The crossover between 
the ballistic and the diffusive regimes depends, in general, on the features of the disorder. However, when the typical 
deflection angle for a single scattering is small, and therefore the transport mean free path £ tr is much larger than 
the mean free path £, a third regime emerges. This regime, known as the directed waves regime, is realized when the 
j — sample width is much smaller than the transport mean free path while it is much larger than the elastic mean free 
path, £ tr >Z>/. In this case, the rays experience many small angle scattering events which result in a diffusive 
dynamics of the ray direction. The total change in propagation direction, however, remains small. 

The focus of our study is on directed waves which are important for many applications ranging from laser communi- 
cations in atmosphere to propagation of acoustic or electromagnetic waves through biological tissues. Similarly to the 
ballistic and the diffusive regimes, the directed waves regime has also been studied in many papers (see for example 
Refs. 5.6.7|H and references therein). However, our results, in many respects, differ substantially from those obtained 
in previous studies. One of the main differences is the slow power law decay of the intensity correlation function in 
space, and the change of its sign, see Fig. [1] This difference affects interpretation of any wave intensity measurement 
which uses a finite aperture apparatus. 

In this article we develop a general method for calculating speckle correlations over distances larger than the light 
wavelength, A. This method, which is similar (but not identical) to the Langevin scheme for the description of classical 
fluctuations^ ' 10 ! 11 , enables one to treat both the diffusive and the directed wave regimes on equal footing. We apply 
the method to the case of directed waves to evaluate speckle correlations and their sensitivity to various perturbations, 
such as a change in the frequency of the wave, a variation of the incidence angle, or a change of the refraction index. 
A short version of these results was published in Ref. [13 . 

The paper is organized as follows. In section Hi Al we present the general method describing speckle statistics. In 
sections III Bl and III CI we consider its limiting cases for angular and spatial diffusion. The treatment of sensitivity of 
speckle patterns to changes in external parameters is presented in section [H Dl In section HTT1 we apply our formalism 
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FIG. 1: The asymptotic behavior of the intensity correlation function, C(p), in the directed waves regime, p is the distance 
between the observation points, A is the light wavelength, £ is the elastic mean free path, Z is the slab width, and 80 and 8 are 
the typical scattering angle of a ray traveling a distance t and Z respectively. 



to study speckle correlations in the directed wave regime, and spatial diffusion. Finally, in section HVl we present our 
conclusions. The derivation of the formalism is deferred to the Appendices. 



II. METHODS OF DESCRIPTION OF SPECKLE STATISTICS 



A paradigm model for propagation of coherent waves through disordered media is the stationary wave equation for 
a scalar field ip(r), 

fcV(r)V>(r) + V 2 V(r) = 0, (1) 

where k = 2n/\ is the wave number, and n(r) = 1 + Sn(r) is the index of refraction. For simplicity we assume Sn(r) 
to be a random Gaussian quantity characterized by zero average, and isotropic correlation function 

(Sn(r)Sn{r')) =g(\r-r'\). (2) 

Here the angular brackets (. . .) denote averaging over the random realizations of n(r). We assume that the isotropic 
function g(r) is characterized by a single correlation length, £ = [J d 3 rr 2 g(r)/3 J d 3 rg(r)] 1 ^ 2 . 

The above model is studied below. The central object of our approach is the ray distribution function, 

/ "'"> = /^/*'* ('"I)** H) < 3 » 

which may be viewed as the density of rays at the point r and time t propagating in the direction specified by the 
unit vector s. In particular the intensity of the wave at the point r is J(r) = |^(r)| 2 = f d 2 s/(r,s). 

The ray distribution function /(r,s) is a random, sample specific quantity whose statistics can be characterized by 
its moments. We focus on the first (/(r,s)) and second (/(r, s)/(r', s')) moments of this quantity. These moments 
quantify the main features of speckle patterns. 



A. General approach to speckle statistics 

In this subsection we discuss a general approach to describe speckles of coherent waves that is valid both in the 
ballistic and diffusive regimes, and holds for a general angular dependence of the scattering amplitude at a single 
scatterer. 

A general method for calculating moments of the ray distribution function is the disorder diagram technique^. If 
I ^> A, and on the length scale |r — r'| > A, this formalism can be reduced to a set of equations for the average distri- 
bution function, (/(r, s)), and the correlation function of the ray distribution function fluctuations, (Sf(r, s)Sf(r', s')), 
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where Sf — f — (/). These equations describe speckles on various length scales: From the ballistic regime to the dif- 
fusive limit, and are similar, but not identical, to the Boltzmann-Langevin equations in the kinetic theory of classical 
particle s 10 ! 11 ' 12 . Thus (/(r,s)) and (5/(r, s)£/(r', s')) can be deduced from the following set of equations: 



a _ 9</(r,s)) = [(/(rj g))] ^ , d2s , w{s _ g/) g , }) _ (/( ^ 



9r 



g . ^/(r,s) _ Jat{g/(r ;S ) } ^£( r , s)i 



(4) 



(5) 



where the integral over the ray directions, s, is normalized to unity, J d 2 s = 1, and the Langevin sources, £(r, s), have 
zero mean and correlations of the form: 



2-7T 

(£(r,s)£(r',s')) = -^(r-r') 



5(s-s')(/(r,s)) / d 2 Siy(s-s)(/(r,s))-(/(r,s))W(s-s')(/(r,s')) 



(6) 



Here W(s — s') is the probability, per unit length, for scattering between propagation directions s and s'. The mean 
free path £ and the transport mean free path £ tr are expressed in terms of W(s — s') as 



r 1 = 



J ds'W{s - s'), t£ = J ds'(l - s ■ s')W(s - s'). 



(7) 



In the Born approximation the scattering probability can be expressed in terms of the refraction index correlator, 
Eq. ©, as 



W(s) 



k 4 



d 3 rg(r)e l 



(8) 



The derivation of Eqs. J4][6]), using the standard impurity diagram technique^, is presented in Appendix lAl On 
spatial scales larger than £, and £t r it is possible to simplify Eqs. ([3][5]) reducing them to a diffusion-type equations. 
Another simplification occurs if the scattering angle at a single impurity is small. Then at lengths greater than the 
mean free path the change of direction of the wave propagation is described by diffusion in the angular space. The 
simplified form of the general formalism in these two limits is considered in Sections IIIBI and III CI 

Qualitatively the form of the correlation function of the random sources, Eq. ([S]), can be understood as follows. 
Inside the random medium the propagating wave can be viewed as a random superposition of plane waves arriving 
from different directions. The relative phases of the different plane waves are uncorrelated. Let us consider scattering 
of this incident wave at a given impurity. Denoting the amplitude of the wave incident in the direction s by i(s) we 
can express the angular dependence of the the outgoing wave, o(s), as 

(.)=i(. )+ a*/^(., s >C), 

where F(s, s') is the scattering amplitude. The intensity of the outgoing wave in the direction s is 

|o(s)| 2 = \i{s)\ 2 -4k J ds'lm [F(a, s')i* (s)i(s')] + 4fc 2 J ds'F{s, s')i(s') . (9) 

The flux into direction s due to scattering, j(s) = |o(s) | 2 — |i(s)| 2 , is a random quantity. Since the amplitudes i(s) of 
the incident wave are uncorrelated for different directions, (i(s)i* (a')) ~ S(s — s')(/(s)), the average flux is given by 



(j(a)) =-4k(f( S ))lm[F(a,s)]+4k 2 J ds'\F(s, s')\ 2 (f(s')) = 4fc 2 J ds'\F(s, s')\ 2 [(f(s')) - (f(s))\, (10) 

in agreement with Eq. ((4]). The last equality in Eq. (|10[) follows from the optical theorem, lm[F(s, s)] = 
kfds'\F{s,s')\ 2 . 

For a specific realization of the incident wave, the flux scattered in direction s differs from its average. In the 
spirit of the Boltzmann-Langevin approach one has to evaluate the fluctuations of microscopic fluxes in the s space 
and substitute them into the kinetic equation as random sources C(s) ~ j(s), see Eq. ([5]). Thus, for the correlation 
function of these quantities, (C(s)C(s')) ~ (j(s)j(s')), and using Eq. ((9]) one gets the estimate, 



(£(s)£(s>))~6(s-s')(f(s)) I d§\F(s,s)\ 2 (f(s)) - (f(s))(f(s>))\F(s,s')\ 2 



in agreement with Eq. ©. Here we took into account the fact that in the limit £ ^> A, F(s, s') the main contribution 
to the flux correlations comes from the middle term in the right hand side of Eq. @ . 
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1. A comparison between Eqs. RJU) and the Langevin description of classical fluctuations. 



It is instructive to compare the method describing classical kinetics of particle o 10 ' 11 ! 12 with the description of 
coherent wave propagating through a disordered media expressed by Eqs. Consider noninteracting particles 

propagating in a scattering medium, and let /(r, s,t) denote their distribution function in phase space. The scat- 
tering process of the particle is random in time and space. This randomness leads to temporal fluctuations of the 
distribution function / even when the incident particle flux is stationary. It is, therefore, natural to decompose the 
distribution function /(r,s;i) into a sum of its average, ((/(r, s; t)}}, and fluctuating part, Sf(r, s;i), characterized 
by the correlation function ((SfSf)). Here ((• • • )) denotes the averaging over time, or over the statistical ensemble^. 

When the elastic mean free path is much larger than the disorder correlation length, £ 3> £, The average distribution 
function satisfies the Boltzmann kinetic equation: 

cot Or 

ds'[/ s ( + ) (s,s') + /ir ) (s,s')]= J A'ff( 8 -s')(((/(r, S ';t)))-((/(r,M)))). (11) 

where c is the particle velocity, I^\s, s') denotes the particle flux from s' to s due to collisions, and /fj (s, s') denotes 
the particle flux from s to s', 

4 +) (s )S ') = W(s-s)/(r, S ';i)(l±/(r lS ;t)), 
lir ) (s,s') = -W r (s-s / )/(r,8,t)(l±/(r,s / ;t)). 

The ± sings in front of /(r, s'; t) correspond to boson (+) and fermion (— ) statistics. Notice however that the quadratic 
terms in ((/(r, s, i))) cancel out in the Boltzmann equation (fTTj) and regardless of the particle statistics. 

The statistical behavior of the fluctuations of the distribution function, 8f(r, s; t), may be deduced from the Langevin 
equation^ ., 

dSf(r,s,t) dSf(r,s,t) - ~ ( 

^ +s ^ =I st {Sf(r,s)} + I L {r,i), (12) 

where It, represents a random Langevin source with vanishing expectation value and two point correlation function 
given by 

{{I L (r, s, t)I L (r', s', f)» - S(t-t')S(r-r') |<5(s - *') J d S "[I { +\s, s") + WO] - [4 + W) + (s, s')] } • 

(13) 

The classical limit of this equation corresponds to ((/(r, s,t))) <C 1. In this case particle statistics are irrelevant. 
The description of the evolution of the average ray distribution function, by the Boltzmann kinetic equation of a 
classical particle holds as long as I ^> £, A. The above formulaes have the following interpretatio n 10 ! 11 : The scattering 
processes are instantaneous and local therefore the correlation function of Langevin sources (|13p is proportional to 
S(r — r')5(t — t'). Thus scattering events generate correlations of Langevin sources that are nonlocal only in the 
space of the particle directions. These are described by the four terms in the curly brackets. The first two terms, 
proportional to 8{s — s'), describe self-correlation generated by flux of particles which scatter from the state s to an 
arbitrary state s" or vice versa. The two other terms in the curly brackets correspond to scattering events from s to 
s', or back. 

The set of Eqs. (I1HI13[) describing the kinetics of classical particle and that of Eqs. (d][S]) describing coherent waves 
have a similar form. We would like to point out significant differences originating from the different nature of fluctua- 
tions. A stationary coherent wave propagating through a disordered sample experiences no temporal fluctuations. In 
this case the spatial fluctuations of /(r, s) result from the random nature of the interference processes associated with 
different quasiclassical wave propagation paths. As a result the random sources, Eq. ©, are (5-correlated in space 
and do not depend on time. In contrast, in the case of classical particles / fluctuates both in space and in time, and 
consequently the random classical sources, Eqs. (I12H13[) are (5-correlated both in space and in time. 

Another significant difference manifests itself in dramatically different sensitivities of these two phenomena to small 
changes of parameters, such as particle's velocities (or wavelength), frequencies, and configuration of the scattering 
potential. In the case of classical particles the correlators ((/)) and ((SfSf)) are insensitive to these changes as long 
as the scattering probability W(s — s') does not depend on the wave length or the energy of the particles. In contrast, 
the coherent speckles exhibit very strong sensitivity to these changes. As a result the form of the correlation functions 
of the random sources describing these sensitivities, see Eq. (|25p . is very different from that in Eq. ITj 
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B. Angular diffusion 

As mentioned above the solutions of Eqs. (|4j[6j) provide description of (/(r, s)) and 5(f(r, s)8f(r', s')) on the resolu- 
tion where |r — r'| > A. A simplified description is obtained when the required resolution is over larger length scales. 
Consider the case |r — r'| ^> £9 , where 8 a = A/£ is the typical scattering angle over a distance of the order of the 
mean free path (notice that the Born approximation implies that 19$ 3> A). The reduction of Eqs. (|4][6]), for this case, 
is similar in spirit to the standard way by which Boltzmann equation is reduced to the diffusion equation. It follows 
from the assumption that /(s, r) changes slowly as function of s on the scale of order Oq. The resulting formulae, 
provided below, describe diffusive spreading of the rays in the space of directions, s. Equation reduces to 

s .^ 1 f> =£)flV 2 </(r>s)>) (M) 



where 



De = \itr (15) 
is the diffusion constant in the space of angles, s, and 

39 sin(0) dtp 

is the gradient operator, with the unit vectors (/)—(— sin <fi, cos <j>, 0) and 9 — (cos <j> cos 9, — sin <j> cos 9, — sin </>) (here 9 
and <j> are the angles associated with polar coordinates). 

The fluctuations of the ray distribution function in this case are described by the equation 

s- dSf !f ,S) = V s [^V s 5/(r,s)-j L (r,s)], (17) 
or 1 

where the Langevin current sources, j L (r, s), are correlated as 



(18) 



(^(r,s)^(f,s)) = J^6 af3 5(s-s)6(r-r), (19) 

where the indices a and (3 denote the vector components in the two-dimensional space of directions that is tangential 
to the unit sphere \s\ = 1. 



C. Diffusion in real space 

If one is concerned with an even cruder resolution, where |r — r'| ^> £ tr , the effective description of the system 
employs the diffusion equation in real space. In this case /(r, s) is assumed to be a nearly isotropic function of s, and 
a slow function of r. Then Eqs. can be reduced to the following set of Diffusion-Langevin equations 2,4 . Namely, 
expressing the wave intensity, J(r), at point r as I(r) — J d 2 sf(r, s), one can reduce Eq. ((4]) to the Laplace equation, 

V 2 (/(r)}=0, (20) 

while the correlator of the intensity fluctuations, SI = I — (I) , can be deduced from the flux conservation condition, 

V-5J = 0, (21) 

with 

S3 = -DVSI + 3 L . (22) 

Here D = £t r /3 is the diffusion constant in real space (notice that according to our convention the diffusion constant 
has dimensions of length). The Langevin current sources, 3 L , have a vanishing expectation value and are characterized 
by the correlation function: 

(J Q L (r)J^(r')) = ^-{I{v))H af} 5{v-v'). (23) 

The boundary conditions for these equations are the conventional conditions for the diffusion equation: SI — at 
open boundaries, and J • n = 0, with n being the normal to the boundary, at closed boundaries. 
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D. Sensitivities of speckles to changes of external parameters. 

The interfering waves travel along different paths, and the lengths of these paths are much longer than the wave 
length. Therefore the phases accumulated along each path are very sensitive to changes of external parameters such 
as the wave number k, the incidence angle of the incoming wave, or a smooth change in the refractive index An(r). 
We will characterize these changes by the control parameter 7(1") = Afc + fcAn(r) where Afc denotes a change in the 
wave number, k. The formalism presented above may be straightforwardly generalized to calculate the sensitivity of 
the speckle pattern to various external perturbations. The sensitivity of the speckle pattern can be characterized by 
the correlator of the ray distribution functions at different values of the control parameter, (Sf(r, s, 0)<$/(r, s, 7)). In 
order to evaluate it Eq. ([5]) should be replaced by two equations. One for df(r,s, 0), and another for 6f(r, s, 7). The 
form of these equations is precisely that of Eq. |5]), however the Langevin sources now depend on the perturbation 
parameter 7. Namely 

s • dSf{r d * ;l) - I st {Sf(r, s; 7)} = £(r, s; 7), (24) 

where £(r, s;7) denotes the Langevin source associated with the value 7 of the perturbation. The average of the 
Langevin sources vanishes. Their correlation function, at different points in space and different values of the control 
parameter, is given by 



(£(r,s;0)£(r',s';7)) = -^(5(r-r')^ 6(s-s')f u (v, s) J d 2 Sl iy(s- Sl )/_,(r, Sl ) - /„(r, s)W(s-s')/-„(r, s') 



v=± 

where /±(r, s) satisfies the equation, 

<9/±(r,s) 



> (25) 



-/ st {/±(r,s)}=±i7/±(r,s). (26) 



dr 

At free boundaries, the boundary conditions for the functions /±(r, s) coincide with the standard boundary conditions 
for the Boltzmann equation. At the boundary with an incident radiation, denoted by S, the functions /±(r, s) are 
determined by the parametric correlations in the incident wave, i.e. 



/+( r > s ) I res = fy^T Jdr'^ {^'i^j i>o ( r +7T 



(27) 



Here the subscript of the wave amplitude "0 denotes the value of the parameter 7. The corresponding equation for 
/_(r, s) is obtained from Eq. [)27p by interchanging the subscripts: 7 <-> 0. 

When the external perturbation is associated with a change in the incidence angle of the incoming wave, Eq. (j25|) 
still holds, however, both /+(r, s) and /_(r, s) satisfy the same equation (fl4"|) . The difference between / + (r, s) and 
/_(r, s) arise from the boundary conditions. (|27l) . We shall elaborate on this issue in Section UlI A 21 

The above formulae describe the speckle sensitivity on the resolution scale larger than the wavelength. As discussed 
in the previous section the formalism simplifies for lower resolution. We conclude this section by providing the relevant 
formulas for the case of angular diffusion, and diffusion is real space. 

1. Sensitivity in the case of angle diffusion 

If the typical scattering angle at a single impurity is small and wave propagation length exceeds the mean free path, 
£, the equation for the fluctuations in the ray distribution function is 



s 



<W(r,s;7) 



= V s [ J D e V^/(r,s; 7 )-j L (r,s;7)] , (28) 



Or 

where the Langevin current sources, j L (r, s; 7) depend on the perturbation 7. These have zero mean and correlation 
function given by 

(#(r,s;0)#(r,S5 7)) - ^'M^U*, «0 s aP S{B-S)S(r-i), (29) 



where f±(r, s) satisfy the equation 

' )J ,r ' S) -D e V 2 J±(r, S )±i 7 f±(r, S ). (30) 



Or 
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2. Sensitivity in the case of the real space diffusion 



Finally, on spatial scale larger than the transport mean free path £ t r, the sensitivity of the speckle pattern may be 
described by the current conservation condition, 

V • SJ = V • (-DV5I + 3 L ) = 0, (31) 

where the Langevin current sources, at different values of the perturbation parameter, 7, are correlated as 

<J Q L (r;0)J^(r'; 7 )) = ^J + (r)I_(r)* a/3 *(r - r'), (32) 

and I± (r) satisfies the equation 

DV 2 I± (r) ± I± (r) = 0. (33) 



III. EVALUATION OF SPECKLE CORRELATION FUNCTIONS AND SPECKLE SENSITIVITIES TO 

CHANGES OF EXTERNAL PARAMETERS 

In this section we shall illustrate the use of the formalism developed in the previous section. To this end we will 
consider the correlation function of speckles and their sensitivity to various perturbations in the regimes of directed 
waves as well as for diffusion in real space. 



A. Speckles in the regime of directed waves 



Consider a situation in which a wave of intensity Iq is incident on a disordered slab of thickness Z, as shown in the 
inset of Fig. [T] The slab thickness is assumed to be much smaller than the transport mean free path and much larger 
than the elastic mean free path, £ tr ^> Z ^> £. Thus rays diffuse in angle, but their total change of direction is small. 
In this regime of directed waves it will be convenient to choose the coordinate system r = (p, z) where z is the the 
direction of the wave propagation in the absence of disorder (8n(r) = 0), and p denotes a two dimensional vector in 
the plane perpendicular to the z-axis. Similarly we decompose the vector of the ray direction as s = (sj_), s z , where 
s z denotes the component in the z direction, while sj_ is a two dimensional vector in the perpendicular plane. The 
rays of directed waves are almost parallel to the z axis and therefore s z « 1, i.e. s sa (sj_, 1). If we denote by 9 the 
typical ray angle at z = Z, then the latter approximation holds as long as 9 <C 1. The results which we present below 
are calculated to leading order in the small parameter 9. 

It is instructive to start with understanding the classical evolution of the average ray distribution function in the 
regime of directed waves. For this purpose we solve Eq. (fT4| for the case where a single ray moving in the z direction, 
impinges upon the slab at the origin r = 0. The assumption that s ~ (sj_, 1) allows one to reduce Eq. (jT5J) to 



d(f(r,s)) 
dz 



<9(/(r,s)) n 9 2 (/(r,s)) 



dp 



- D, 



ds 2 ± 



= 0. 



The boundary conditions are 



(f(r,s))\ z=0 = i 5(p)6(s 



where the amplitude io denotes the incident ray intensity. The solution of the above problem takes the form 



</(r,s)> = 



3i 



47T 2 L>2 Z 4 



exp 



3p 2 3s_i_p 



D g z A D e z x 



Dgz 



(34) 



(35) 



(36) 



Dgz. It also shows 



It demonstrates the diffusive behavior of the ray direction as it propagates in the slab, 
that deviations in real space grow in a superdiffusive manner 19 , p 2 ~ Dqz 3 . 

After this preliminary consideration we turn to study intensity correlations of directed waves. To be specific we 
consider a plane wave (not restricted by a finite aperture) incident on the disordered slab in the z-direction. In this 
case the average ray distribution function is independent of the perpendicular coordinate p and can be easily obtained 
by integrating Eq. (|3"6"|) over p, 



(/M = 



^0 

AnDez 



exp 



4D g z 



(37) 
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The intensity correlation function, 



C(cfr) = (SI{r)SI{r + Sr)), 



(38) 



where 51 (r — J(r)— (/(r), is independent of the transverse coordinate and depends only on the propagation distance Z 
and the difference coordinate Sr. The behavior of this correlator as a function of Sr = (p, Sz) is strongly anisotropic. 
Consider first the case where the observation points are located along the z axis (i.e. p = 0) near the point z = Z. 
In this case we obtain 



C[Sz) = 



T 2 



4k 2 9 i Sz 2 ' 



(39) 



where 9 = \JDqZ is the accumulated scattering angle, and the condition I <C Sz <C Z is assumed. This formula, 
which also approximates the behavior for nonzero p as long as Sz p/9, matches the results for the diffusive case^il, 
Z 3> Itr when 9 is of order unity. 

A more complex behavior of the correlation function appears when Sz < p/9, i.e. when the observation points are 
located essentially in the plane perpendicular to the z-axis. A general formula for C(p), in this case, is derived in 
Appendix B. The expression takes the form 



C(p) 



i 



ADgk 2 



Z-Jt 



c 



(40) 



where g(p) = J dzg(^J p 2 + z 2 )/ J dzg(z), and Jq(x) is the Bessel function of zeroth order. 

The integral in Eq. (|40[) contains a term proportional to a (5-function, 2D £ 5z 5(p). This term represents the rapidly 
decaying (at p ~ A/0) part of the correlator. It results from the scmiclassical approximation employed in the derivation 
of Eqs. (J4J6]), which limits the spacial resolution to p > A. In order to resolve the spatial structure on smaller scales 
some of the diagrams discussed in Appendix [A] should be calculated more accurately. The result of this calculation 
shows that the S function contribution to the correlator C(p) is in fact a contribution of the form J ( 2 e~ 2 ( feep ) where 
9 2 =D e Z. 

As we show below, C(p) contains also a slowly decaying term. The latter, which has been overlooked in previous 
studies, clearly has important consequences. In order to understand this term it will be instructive to explain, first, 
the origin of the short ranged contribution to C(p). As we show now, it arises from a superposition of statistically 
independent contributions of waves moving in all possible directions. Let us assume that wave function at a given 
point on the screen is a sum of plane waves. The distribution of directions of these plane waves is dictated by the 
diffusive nature of the rays in the system, thus 



(41) 



where s_l „ denotes the direction of the z^-th contribution and A v is the corresponding amplitude. We shall assume 
that A v are statistically independent variables, with zero mean and fluctuation strength given by 



(l^l 2 ) = 



4nD e Z 



(42) 



The average (lA^I 2 ) may be interpreted as the "classical" probability to find a plain wave moving in direction sj_ v . 
It may be obtained from the solution of Eq. (|34[) with boundary conditions which correspond to an impinging plane 
wave of density 7 , (/(r,s))| z=0 = Iq5{b±). 

The above assumptions imply that, at a given point in space, ip(p) is approximately a Gaussian random variable, 
as a result of the central limit theorem. Moreover, the wave function at two different points, tp(p) and i/j(p'), are also 
described by a joint Gaussian distribution function provided the distance between these points is sufficiently small 
such that one may assume that the same set of wavelets arrive to both points. 

Assuming the observation points p and p' to be sufficiently close to each other, consider the ensemble average 
{I{p)I{p')) = (ip(p)ip*(p)ip(p , )ip*(p'))- Using the fact that within a small vicinity of space, ip*(p) may be considered 
as a random Gaussian function, one deduces that {ip{p)ij}* (p)) (ip(p')tp* (p')) + (t/j(p)t/j*(p'))(ip(p')i(}*(p)), and hence 
the density correlation function is given by 



c{p-p') = m P w{p'))\ 2 



(43) 
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Now from (|41 [) and the statistical independence of the amplitudes A v we see that 

(v>(p)v*(p')> = (y2A v At,e ik ( e ^p- s ^'P'^\ =Y,{\M 2 v ks ^ p ~ p,) ■ (44) 

\v,v' I V 

The replacement of the above double sum by a sum over one index is equivalent to the assumption that the interference 
terms of different amplitudes average out to zero. This traditional procedure in semiclassical analysis, known as the 
"diagonal approximation" leaves only the classical contribution. Thus substituting (|42p in (|44|) and replacing the sum 
over v by an integral over sj_ we obtain an expression for (ip(p) , ip*(p')}, and from (|43[) we conclude that 

C{p-p')^lte- 2 ^P-p'^. (45) 

This expression, which shows very fast decay of correlations on a scale of order X/9, has rather limited range 
of applicability. The reason is that the description of the wave function in the sample, using the superposition of 
independent plane waves (|4ip . gives reasonable approximation only when the observation points are very close. At 
larger distances diffraction and quantum impurity scattering give rise to correlation of rays which manifest themselves 
in a slow decay of the intensity correlations, as well a change of sign. These effects are described by Eq. (4Ql and 
illustrated in Fig. 1. At various spatial separations p one can obtain the following asymptotic expressions for the 
intensity correlator: 

e' 2{ - k0 p) 2 if p ~ aX/9, 

WWmTp if aX / d « P « i0 o, 

b - D e /3 ;f oa.^ n ^a 7 (46) 



CM 
I 2 o 



if ee < p < ez 

b 3P 2 e --£f?? if H7. n ZfL 

8o 



where a 2 — \og(k£9 3 / 9q) , b\ = dxg(x)/(, is a constant of order unity, bi — 3 1 / 3 r(5/3)/8 w 0.163, and 63 = 
27/128 w 0.21. The qualitative form of the function C(p) is shown in Fig. [1] 

In order to clarify the connection between ray diffraction and the slow decay of the density correlations, let us 
focus on the regime £9q «/)« 0Z . Consider two points separated by a distance p. The correlations of the wave 
intensity in these points emerge from coherent waves which simultaneously arrive to the two points. These can be 
generated by diffraction which acts as a beam splitter, and modeled by the Langevin sources in Eq. ((5|). Now, the 
superdiffusion nature of the ray dynamics in the sample implies that the relevant points where diffraction takes place 
should be located at distance of order Az from the screen, where p 2 — DgAz 3 . The wave intensity emitted from these 
diffraction points decays as 1/Az 2 , and therefore the correlations which they generate are proportional to Z?^ 3 /p 4 / 3 . 

The above crude argument explains the power law decay of C(p), in the regime £9q <C p <C 0Z. Yet a closer 
examination of the integrals leading to this results shows that the contribution from diffraction points (or Langevin 
sources) that are closer to the screen than, Az — (Dgp 2 ) 1 / 3 generate anti-correlations, while those that are at larger 
distances provide positive correlations. This behavior may be expected since diffraction points located too close to 
the screen generate rays which may arrive to either one of the observation points but not to both of them, therefore 
they lead to anti-correlated behavior. On the other hand, coherent waves generated by diffraction that took place at 
distances larger than Az, get, in general, to both points, and therefore generate positive correlations. 

From this picture, and the finite width of the slab, it follows that for sufficiently large distance between the 
observation points, p 2 ^> DgZ 3 = (Z9) 2 , diffraction events can generate only anticorrelations. Thus C{p) must 
experience a sign change in the vicinity of p = 9Z . 

Finally, we mention that the tail of the correlation function (the regime p > Z9 2 /0q) is also described by Eq. (|40[) . 
However, it depends on the precise form of the disorder correlation g(r), since this limit is dominated by rare scattering 
events. 

The power law nature of the density correlations of directed waves have important consequences regarding the 
statistics of the signal measured by sensors with large apertures compared to the wave length. Let 

P = J d 2 pd 2 s n- s/(r,s), (47) 

denote the signal measured by the sensor, where n is a unit vector perpendicular to the sensor surface, and p is a two 
dimensional vector which parameterizes the sensor surface. If the sensor aperture is circular, with radius R, and its 
surface is perpendicular to the propagation direction, i.e. s ■ n ~ 1, then the integrated power measured by the sensor 
may be approximated by an integral over the wave density 

P(R) = I d 2 pl(r), (48) 
J\p\<R 
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where r = (p,z). The random fluctuations of J(r) imply that P(R) is also a random quantity. Its average may be 
expressed as an integral over (I(r)), while the variance of its fluctuations is given by 

((SP(R)f) = / d 2 pd 2 p'C{v-v'), (49) 

J\p\,\p'\<R 

where C(r — r') is the density correlation function (j4"0"|) . 

Clearly, the fluctuations of P{R) strongly depend on the slow power law tails of the correlation function as well 
as its sign change. The asymptotic behavior of the variance of these fluctuations, for a circular sensor with aperture 
radius R is given by 

! 7r _r_ b 1 R g\ T? ^ Pfi 

2k 2 8 2 k 2 6 i t8 ' g <^ri<^lU , 

„ i b' 2 (D eR )^ pa^u^OV (50) 



where &' x = 26itt/3, b' 2 = 3 4 / 3 r(5/6)7r/2 n / 3 r(7/6), and b' 3 = ^/Jj2^. 



1. Speckle sensitivity to change of the wave frequency 

Consider the sensitivity of the speckle patterns of directed waves to a change in the wave frequency: Auj = cAfc, 
where c is the speed of the wave, and A; is the wave number. Using Eqs. (124II26|) . with the appropriate control 
parameter, 7 = Ak, treated on a perturbative level, one may identify the scale of the change in the control parameter, 
where the new speckle pattern essentially lost its correlations with the initial one (i.e. the speckle pattern at 7 = 0). 
For the wave frequency perturbation this scale is found to be 

= ^ (51) 

A qualitative explanation of the scale ui* is similar to that given for the sensitivity of the conductance 
fluctuation o 15 i 20 . Let us estimate the characteristic change in the phase of a typical orbit due to the frequency 
change Aui: The typical length spread of the orbits, in the directed waves regime, as follows from their superdiffusive 
nature, is of order 2 Z. Therefore the phase difference between a given orbit and the same orbit different frequency, is 
of order AkZ9 2 where Ak = Aw/c is the change in the wave number. Thus a complete change of the speckle pattern 
occurs when the phase, AujZQ 2 jc is of order one, namely Auj ~ c/ZO 2 ~ co* , in agreement with Ref. [8|. 



2. Sensitivity of speckles to change of the angle of incidence 

Consider the case where rays propagate through a disordered slab whose one edge is located at z = 0. A plane 
wave, moving in direction approximately parallel to the z axis, impinges the slab, at z = 0. The speckle pattern 
formed on the second edge of the slab, at z — Z, will be sensitive to the precise angle, <f>, of the incoming wave. The 
latter takes the form yj = t/I^ exp[ikz cos <f> + ikpsm.<f>]. 

As mentioned in the previous section, the sensitivity in this case is characterized by the correlation function (|25|) 
of the Langevin sources (C(r, s; 0)£(r', s'; <fi)}, where both /+(r, s) and /_(r, s) satisfy the same equation 

s d/j^s) 
or 

However their boundary conditions are different. They are determined by the Wigner transforms of a product of the 
incoming wave parallel to the z axis, by the complex conjugate of an incoming wave at angle ±<f> (evaluated at z = 0). 
Thus the boundary conditions for Eq. ([52"]) are 

f ± (p, z = Q) = I e ±lks ^S(s - s ), (53) 

where Sq — (cos0, sj_) w (l,sx) denotes a unit vector in the direction of the incoming wave, and |sj_| = sin</> 0, 
assuming <f> <C 1 • Solving the above equations one can identify the characteristic scale for the change in the incidence 
angle: 



^ = kZB- 



(54) 
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This result has simple interpretation. Consider a given point on the screen. The wave intensity at this point is 
determined by the interference of all the rays which originate at z — and reach the same point. The nature of the 
ray dynamics, in the directed wave regime, implies that the the original distance between two rays which reach the 
same point at the screen is of order of ZO. Now if we change the incidence angle by some small amount <fr <C 1, the 
phase difference between two such rays is of order kZ8(f>, where k is the wave number. The interference of these rays 
will be completely different when this phase difference is of order one, i.e. kZ8(j>* ~ 1. From here we obtain (fS"4|) . 



B. Speckle statistics in the diffusive regime, Z ^> t tr 



In what follows we complete the picture of speckle statistics by presenting the well known results of speckle cor- 
relation functions and sensitivities for the diffusive regime, Z > 4. For simplicity we consider the situation where 
I = £tr, and set the resolution scale to be larger than the wavelength, A. Furthermore, as in the previous section, 
we shall consider the infinite slab geometry shown in Fig. [TJ and assume that plane wave, moving in the z direction, 
impinges the system at z = 0. 



1. The intensity correlation function 

Our first step is to solve Eq. (|20|1 for the average intensity. The boundary conditions in this case are I(r)\ z —z = 
and dI(r)/dz\ z= o = —Jq/Z, where Jo is the flux of the incoming wave, and D is the diffusion constant. Thus 

(/(*)) = J £^£. (55) 

This solution implies that the flux inside the sample is (J z ) = —Dd(I(z)) / dz = Jo£ tr /Z and therefore the average 
transmission coefficient through the slab is ratio of the mean free path to the width of the slab: 

(T) = = L (56) 
Jo & 

Notice that in our conventions the diffusion constant, D has dimensions of length, and is the transport mean free 
path. 

Consider now the density correlation function (|38|) . Solving Eqs. (|2Tj) . (|22[) . and calculating C(r), using the corre- 
lation function of the Langevin sources (f!?5)) (evaluated with the help of (|55|) ), we obtain 

f A<r<£ tr 
C(r) = {I(z)f I , (57) 

{ im^r ^r^Z 

where it is assumed that the observation points are far form the end of the slab, i.e. < ir « z « Z. 

The above result shows a power law decay of the speckle correlations which is similar to the case of directed waves. 
Yet, unlike directed wave, the transmission coefficient of the system in diffusive systems experience sample specific 
fluctuations. This is due to the finite amount of backscattering which can be safely neglected in the case of directed 
waves. In order to evaluate the magnitude of these fluctuations, let us consider the integrated flux passing through 
the slab: 



SJz = ^ I d 3 rJ z , (58) 
<v 



z 



where V denotes the volume of the slab. Here we assume the slab to be finite with dimensions X, Y and Z, such that 
Z <C X, Y. Now, as follows from Eq. (f!Z2]) . the current J z contains two contributions: 

J Z = -D^-6I+JL. (59) 
oz 

The first contribution vanishes upon integration over space, therefore the fluctuations in the total current are essentially 
due to the contribution from the Langevin sources: 

<(^) 2 > = i /^ 3 ^V(jf(r)J 2 V)>- (60) 
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Substituting Eq. (|23]) for the correlation function of the Langevin sources, and evaluating the integral we obtain: 

((6J M )>) = V™ (61) 
From here we conclude that the fluctuations in the transmission coefficient scale as: 

(T) 2 (J Z ) 2 V 2 18n£XY N' K ' 

where N ~ vV/tq ~ £XY/X 2 Z is the total number of eigenfrequencies lying within frequency band of width 1/tq, 
centered at the frequency of the incoming beam. Here v ~ 1/cA 2 is the density of states of the slab (per unit volume), 
tq = Z 2 jlc, is the typical time of diffusion through the sample, and c is the wave velocity. 



2. Sensitivities of the speckle pattern in the diffusive regime 

Below we summarize the results of the speckle pattern sensitivities to various perturbations in the diffusive regime. 
These results are obtained by solving Eqs. (|3H33|) and identifying the the relevant scale of the perturbation parameter. 
The sensitivity to a change in the wave frequency is characterized by the frequency scale of the order of 

cl , s 

= -p, (63) 

where c is the wave velocity, and £ is the elastic mean free path. This frequency scale is the inverse time which takes 
the wave to propagate through the sample. 

The sensitivity to a change in the angle of the incoming wave, <fi, is characterized by the scale 



1 



kJlZ 



(64) 



The interpretation of this result is similar to that presented for directed waves. Here, however, the diffusive nature of 
the ray dynamics implies that the the original distance between two rays which reach the same point at the screen is 
of order of \flZ. Therefore the interference of these rays will become completely different when the phase difference, 
due to the change in the incidence angle, is of order one, i.e. ky/tZ<j>* ~ 1. This condition leads to (|64p. 

Finally let us discuss the sensitivity of the transmission coefficient to a change in the angle of incidence, in a finite 
three dimensional system. This sensitivity may be described by the correlation function of the fluctuations ST (8) at 
two different angles, and the result takes the forrn^ 

{5T(e)5T(8>)) f &J±W\ if i<|0-0'l<i 
( ST2 ) IU if \<\0-0'\. 

As we show above, the fluctuations in the transmission coefficient follow from the fluctuations in the current due to 
the Langevin current sources. Therefore, one expects that the above correlation function can be deduced from the 
correlation function of the Langevin sources (|32[) where 7 stands for the change in the incidence angle of the incoming 
wave. This procedure, indeed, gives the result within the range 4 < \9 — 0'\ < 4. However for larger difference 
in the angle of incidence, i.e. j < \8 — 8'\, the behavior is dominated by an additional contribution which is not 
described by the Boltzmann-Langevin approach. This contribution can be calculated from a diagram which contains 
two Hikami boxes, as shown in Fig. [2l In real space it may be associated with pair of orbits which intersect twice 
during their propagation in the system. 

At this point it is instructive to mention the relation between Eq. (|65p and the universal conductance fluctuations 
of mesoscopic metals. The conductance in these systems is proportional to the integral of the transmission coefficient 
over the angle, G oc J T(6)d0. Therefore according to Eq. (l65|) the main contribution to the conductance fluctuations, 

((SG) 2 ) cx J d0dO'5(T(8)ST(O')), (66) 

comes from the interval of large angle difference, j < \9 — 6'\. As a result we have (SG) 2 ) ~ e 4 /h 2 
dimensional system where all dimensions are of the same orde r 16 i 20 . 
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FIG. 2: The diagram contributing to the transmission coefficient correlations at large angles. The gray boxes are Hikami boxes, 
while solid lines connected by dashed lines represent averaged Green functions, see appendix A for details 

IV. CONCLUSIONS 

We have developed a method of description of speckle statistics in elastically scattering media which can be applied to 
both diffusive and to the ballistic regime. Our main result is given by Eqs. (J4J6J) , which have a form of kinetic equations 
with random sources. Though the derivation of these equations in Appendix [Al involved the Born approximation for 
the amplitude of scattering on individual scatterers, we believe that the region of the applicability of these equation 
is much broader: They are valid as long as the Boltzmann kinetic equation Eq. ((4| holds. Namely, I 3> A, £, and 
|r-r'| » \i. 

We would like to mention that the results presented above substantially differ from those known in the literature 
(see for example Refs. 5.6,7.8). First, the correlation function (|40|) exhibits a universal long range power law behavior 
over a wide range of distance, p. The only non-universal regimes are at the tail, p 3> Z9 2 /9 , and the short distance 
region, p ~ £ . This result is in contrast with the results presented in Refs. 0000 where the intensity correlator 
C{p) depends on the detailed form of the disorder correlation function, g(r), and usually decays exponentially at 
p > £. Second, C{p) changes its sign as a function of p. This property is a consequence of the current conservation 
and it is absent from previous studies. For instance, the sign change of C(p) implies, that the fluctuations of the 
integrated intensity over disks of radius R > ZQ is proportional to R, see Eq. (jSTij) . rather than R 2 , as would follow 
from Refs. These differences will affect interpretations of any measurement of speckles done with the help of 

a sensor aperture that is much larger than the wavelength. 

Finally we would like to mention that our results may be easily extended to cases with light polarization, optically 
active media, Faraday effect, and coherent short wave pulses as long as their duration is longer than r = £/c. These 
issues are left for future studies. 

This work has been supported by the Packard Foundation, by the NSF under Contracts No. DMR-0228104, and by 
the Israel Science Foundation (ISF) funded by the Israeli Academy of Science and Humanities, and by the USA-Israel 
Binational Science Foundation (BSF). 

APPENDIX A: DERIVATION OF THE MAIN EQUATIONS 

The derivation of Eqs. ([4]), (O, and ([6]) is based on the standard impurity diagram technique^. And the relevant 
diagram blocks were derived in numerous works. However in most cases the calculations were done either for the 
case of delta-correlated disorder potential or in the diffusive regime. In this paper we deal with a general situation 
of an arbitrary angular dependence of the scattering cross-section. Therefore below we outline the derivation of our 
formalism and present expressions for the main diagram blocks. 

The wave equation (fT]) can be written in the form of a stationary Schrodinger equation for a particle moving in the 
presence of a random impurity potential, 

V(r) = -2k 2 5n(r). (Al) 

The solution of Eq. |T]) can be written as -0( r ) = fdr'G R (r, r') J(r'), where J(r') is the source of radiation and 

G r / a {y 1 y i ) is the retarded Green function, G R I A = (jt 2 + V 2 — V ± ir^j . Here V denotes the impurity potential 

operator. This reduces the problem of speckle statistics of coherent waves to that of averaging products of retarded 
and advanced Green functions. The latter problem can be treated using the impurity diagram technique™ We derive 
the expression for the various diagram blocks below. 
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FIG. 3: The self-energy diagram in the Born approximation. 

1. Average Green function 

In the Born approximation 21 the self-energy, £(fc,p), is given by a single diagram in Fig. [3] 
Its evaluation gives for the disorder-averaged Green function, 

° R/A ^ = k>-/±ikl-i > (A2) 
where the mean free path is given by Eqs. 10) and ([8]). 

2. Derivation of the Boltzmann equation 

To derive the Boltzmann equation we will need to evaluate products of Green functions at two different frequencies 
corresponding to wave numbers, k± = k ± Sk/2. 

The spatial evolution of the ray distribution function, Eq. ([3]), can be obtained by expressing the solution of the 
wave equation in terms of the Green's functions and performing disorder averaging. In the leading approximation in 
\/£ the ray distribution function evolution is described by the sum of ladder diagrams. 

Each disorder-averaged Green function is strongly peaked in the momentum region where the on-shell condition is 
satisfied, k — |p|. In the limit of dilute scatterers the width of this peak, ~ X/£, is much smaller than the typical 
momentum momentum transfer at each collision. Therefore the integration over the magnitude of momenta in the 
Green functions can be carried out separately and before the direction integration. Defining s as the unit vector along 
the momentum p = ps we evaluate the product of the disorder-averaged retarded and advanced Green's functions 
integrated over p, 

^ 4 W ^rTG R (k + ,ps + q/2)G A (k-,ps-q/2) = — —- —. (A3) 



Bsk,q(s) Jo 27r 2 ' ' -i6k + isq + £ 1 ' 

The ray distribution function /(s, q) is given by the sum of ladder diagrams and can be expressed in a compact 
way using the operator notations, 

oo i 

f = (b^w)" fo ={B-w) Bf . (A4) 

n=0 

Here /o is the initial ray distribution function, B is the integral operator whose kernel in the Fourier representation 
is given by Eq. (|A3I) . and W is the integral operator acting in the space of directions, 



Wf(s) = [ ds'Wis - s')f(s'), Wis - a') = — w(k[s - s']). (A5) 

J 47T 

Multiplying Eq. (|A4[) by (^13 — W^j from the left and using Eq. (0 we obtain the Boltzmann-Langevin equation, 

(-iSk + s-V r ) f(s,r) - I st [/] = A (A6) 
C=(-iSk + 8-V r + r 1 )f(s,r), (A7) 

where the collision integral I st [/] is defined in Eq. (j4|) . 

If one is interested in the average ray distribution function then /o in Eq. (|A6|) should be understood as the ray 
distribution function of the incident radiation at the boundary of the disordered medium. In this case the source 
vanishes in the interior of the medium and the average ray distribution function satisfies the usual homogeneous 
Boltzmann equation. 
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FIG. 4: The diagram for the irreducible correlator of ray distribution functions at different points. Two impurity ladders 
emanating at the radiation source enter the Hikami box, represented by the hashed octagon, from left and right. The ladders 
going to the observation points leave the Hikami box from the top and the bottom. 




FIG. 5: The three diagrams for the Hikami box; H w , # (2) , and if (3) . The ladders coming from the radiation sources enter 
the Hikami box from left and right and are characterized by the four-momenta g 3 / 4 = (W3/4, Q3/4) and the ray directions S3/4. 
The ladders going to the observation points exit the Hikami box from the top and the bottom and are characterized by the 
four-momenta q 1 / 2 = (1^1/2,91/2) an d the ray directions S1/2. Note that in our notations the Hikami box contains a single 
impurity line for each of the incoming ladders and no impurity lines for the outgoing ladders. 



3. Hikami box 



Next let us consider the diagram in Fig. 0] that represents the irreducible correlator of the ray distribution functions. 
It allows the following interpretation which is at the heart of the Boltzmann-Langevin approach developed in this 
paper. The impurity ladders connecting the observation points to the Hikami box propagate the fluctuations of 
the distribution function from the Hikami box out to the observation points. This propagation is described by the 
inhomogeneous Boltzmann-Langevin equation (|A6[) . Then the right hand side of Eq. (|A6[) may be interpreted as the 
"Langevin force" that results in the fluctuations of the ray distribution function. The fluctuations of the Langevin 
force are described by the Hikami box connected to the impurity ladders going out to the radiation sources. Since 
the latter define the average ray distribution function we see that by evaluating the Hikami box we will relate the 
fluctuations of the Langevin force to the average ray distribution function. 

The Hikami box with the external legs is given by the three diagrams in Fig. [5J It is characterized by the four- 
momenta qi = (u)i, q^) and the unit vectors Si characterizing the ray directions. Here i = 1, 2 correspond to outgoing 
momenta (ladders going to the observation points) and i = 3, 4 to the incoming ones (ladders coming from the radiation 
source). The momenta satisfy the conservation law, q\ + q 2 = q 3 + q±. The analytic expression that corresponds to 
the first diagram (with no impurity line) is 

#P = &{si-s 2 ){An) 2 W{ Sl ~ s 3 )W(s 2 - Si ) J 2^-G R (k + u; 1 ,p + q 1 )G A {k, P ) 



xG R (k + uj^p + q i )G A {k + uj 1 - uj 3 ,p + q 1 - q 3 ) 



k' 2 



5(si - s 2 



W{ Sl -s 3 )W(s 2 -s 4 ) ( 1 



B1B2 



1 

Bl 



(A8) 



where we used the shorthand notation Bi 
B 1 +B 2 =B 3 +B i . 



Bski,q ( s ) (with s = Si = s 2 ) and utilized the momentum conservation, 
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The second diagram in Fig. [5] contains an impurity line connecting the two advanced Green functions (between q 3 
and q 2 , and q^ and q\ respectively). It is given by the expression, 

^k ({*},{*}) = {^fW{s 1 -s 3 )W(s 2 ~s 4 )W{s 1 -s 2 ) J ^-G R (k + u; 1 ,p + q 1 )G A (k 1 p) 

xG A (k + LU 1 -cu 3 ,p + q 1 -q 3 ) J P^fG A (k,p')G R (k + LU i ,p' + q i ) 



xG A (k + uj 1 -w 3 ,p' +q 1 -q 3 ) 

7T Wjsi - s 3 )W(s 2 - 8 4 )W{8 1 - gg) 
k 2 BiB'JSaP'i 

where unprimed B's depend on s\ and primed ones on S2, 



(A9) 



Bi = B Uuqi {s x ), (A10) 
B'i = B Ui , q .{s 2 ). (All) 

The third diagram of the Hikami box contains an impurity line connecting the two retarded Green functions 
(between qi and q 3 , and q 2 and q± respectively). It is given by 

4 3) ({»}»{«<}) = (4^) 3 W>2 - s 3 )W^! - s 4 - s 2 ) J^.G R (k + iJ 1 ,p + qi )G A (k,p) 

xG R (k + uj il p + q i ) J ?^G R (k + iJ 1 ,p' + qi )G R {k + uo A ,p' + q 4 ) 



xG A {k + wi - ws,p' + q t - q 3 ) 

n W(s 2 - s 3 )W{ Sl - s i )W(s 1 - 8 2 ) 
k 2 BiB 2 B' 3 Bi 



(A12) 



Next we make use of the fact that in Eqs. (|A8[) , (|A9[) . and (|A12|1 the operators with indices 3 and 4 act on impurity 
ladders 3 and 4 that go out to the radiation sources. These ladders are equal to the average ray distribution functions, 
(fa)- In the interior of the medium the latter obey the Boltzmann equation (|A6p with the vanishing right hand side, 
see discussion below Eq. (|A7j) . Therefore we have 

B- x W{f s ) = (fs). (A13) 

Using Eq. (|A13|) and combining Eqs. (|A8[) . (|A9|) . and (|A12|) we obtain the correlator of the Langevin forces that 
enter the right hand side of Eq. (|A6|) . As a result we can describe speckle fluctuations in the framework of the 
Boltzmann-Langevin scheme, Eqs. (|5|) and ©. 



APPENDIX B: DERIVATION OF FORMULA pO)) 

In this appendix we derive formula (|4"tJ)) for the intensity correlation function in the directed waves limit. For 
this purpose we employ the parabolic and the Markov approximations. Namely the scalar wave equation |1| is 
approximated by a simpler equation, obtained by substituting ij} — * e lkz ip(r) into Q]) and neglecting second order 
derivatives of the wave function with respect to z. The resulting equation takes the form of a Schrodinger equation 
where the coordinated associated with the propagation direction, z, plays the role of fictitious time: 

dib 1 / d 2 d 2 \ 

i i = -T k (M + w) i)+kSn{r ^ (B1) 

The analysis of this equation is further simplified when the Markov approximation is employed. The latter corresponds 
to the situation where the disorder correlation function is anisotropic: It is delta correlated in the propagation direction 
z, and long ranged correlated in the perpendicular directions: 

(6n(r)Sn(r')) = g ± (p - p')6(z - z'). (B2) 

Here angular brackets denote disorder averaging, and g±(p) represents the disorder correlation function in the p = 
(x,y) space. We shall assume that 5n(r) is gaussian random function and that g±{p) is isotropic. 
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These approximations however, do not imply, necessarily, diffusive motion, and therefore applies also for length 
scales shorter than the mean free path. Within these approximations the Green function associated with, Eq. Q, 
henceforth called "diffuson" and denoted by T>(p, p; z), satisfies an equation of the form: 



d p d \ f d 2 q 

dz + kdp) P ( p ' p;z ) ~ k J 4^2^ q ) (PiP'P-^z) -^(PiP; 2 )) = $(P- Po) s (p- Po)tf(z), 



(B3) 



where gj_(q) is the Fourier transform of gx{p)- Notice that here the momentum p is a two component vector in the 
space perpendicular to the propagation direction. 

The above equation can be simplified by Fourier transforming it with respect to the momentum, p. Thus if denote 
by x the variable conjugate to the momentum p, and V denotes the Fourier transform of the diffuson V, then 



d id 2 



dz k dpdx 



V-k 2 [g ± (x) -gx(0)]T> 



ipox 



4tt 2 



■Hp- Po)$( z )- 



(B4) 



To solve this equation we further take its Fourier transform with respect to z and p (with conjugate variables denoted 
by q z , and q respectively): 



iq z 



q d 

k dx 



V - k 2 \g ± (x) - gx(0)]T> 



e ip x-iqp 
4^2 ' 



(B5) 



Now let us decompose the vector x into its components: xm parallel to the vector q, and x± perpendicular to that 
vector. Then the solution of the above equation takes the form: 



V 



ke- i( lPo 
4n 2 q 



dx' e iponX '\\ +ipo±x - 



k 

exp — 

q 



dx'{ (iq z - k 2 (gx&^xx) - 9x(P))) 



(B6) 



where under the assumption of isotropy in the plane perpendicular to the propagation direction g±(x\\,x_i 



g ± \ x /x 2 +x 2 



Now, taking the inverse Fourier transform with respect to q z , integrating over a;|| , and Fourier 
transforming the result with respect to x we obtain the result for the diffuson: 



D(p,p;p ,p ;z) = /0/ e ^(P-P„) +J Po(x-t 9 )- i px exp ^(gx(xlxx)- 9 x(0)) 



(B7) 



If we assume boundary conditions where the average distribution function, at z = 0, is given by fo(p,p), then for 
z > the average distribution function is given by the integral: 



f(p,p,z)= / d 2 p d 2 p V(p,p; p ,p Q ; z)f (p ,Po) 



(B8) 

In particular assuming the incident wave, at z = 0, to be a plane wave pointing at the z direction, fo(p, p) = 47r 2 /o<5(p) 
where Iq is the density, the above integral reduces to 



/(p, z) = Jq J d 2 xexp [-ipx - 2k 2 (gx(x) - gx(0))] 



(B9) 



This formula is exact assuming the parabolic and the Markov approximation. Namely it holds as long as / ^> £ 
(Markov approximation), and £ 3> A (small angle scattering, i.e. parabolic approximation). It holds for any distance 
z < Itri and for any value of the momentum p. It may be further simplified if we assume z>f where £ is the elastic 
mean free path. In this case the dynamics is of diffusive nature in the angle of directions and one may approximate 
the correlation function gx(x) using Taylor expansion near x = 0: 



gx{x) ~gx(0)-D e x 2 



(BIO) 



where Dg = — <?" (0)/2 {g'[_{x) denote the second derivative of gx{x) with respect to x) is the angular diffusion constant. 
Substituting (|B10|) into (|B9[) and preforming the integral over x yields: 



nip 

k 2 D 9 Z 



exp 



P 

Ak 2 D 9 z 



z > I 



(Bll) 
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Let us now consider the fluctuations of the distribution function. Using Eqs. (O and ([6]), one may write their 
corresponding correlation function as 



(5f(p,p,Z)Sf(p',p',Z)) = J Z dz J d 2 p"d 2 p 1 d 2 p 2 D{p,p;p",p 1 ;Z - z) 

x ( Sx (0)<J(pi - P2) - g±(pi - P2)) /(Pi, «)/(P2, z)V{pf, p'; p", P2 ; Z 



(B12) 



where, as before, this result has been obtained under the parabolic and the Markov approximations. The density 
correlation function, C(p — p' can be deduced from (|B12j) by integration over p and p': 



C(p-p') = J d 2 pd 2 p'(5f(p,p,Z)Sf(p',p\Z)) 

Thus substituting (|B12[) and (|B7|) into the above formula, a and performing the integral over p" yields 

d 2 q d 2 pi d 2 p 2 



C{p) = / dC 



47T 2 47T 2 47T 2 



x exp 



(ffx(O)J(pi - pa) - g±(pi - P2)) /(Pi, C)/(P2, C) 



(B13) 



(B14) 



iq p 



Pi - P2 



(Z-0 )+2k 



z-c 



dr,(g(r,q/k)-g(0)) 



This formula is obtained essentially by introducing one Hikami box into the diagrams. The small parameter controlling 
this approximation (i.e. the neglect of additional Hikami boxes) is I ^> £, 2 /X. The distance between the observation 
points should be larger than the disorder correlation length, \p — p'\ > £. 

The above integral can be further simplified if we assume the width of the system, Z , to be much larger than the 
elastic mean free path, Z 3> I. In that case, as can be seen from formula (|B11[) . the width of the average distribution 



function /(p, £) at £ ^ £ is much wider than the width of g±(p), as the width of first function is ky/Dg( ~ jyCJI, 
while the second is of order Therefore, assuming the integral over z to be dominated by points near the screen 
z = Z (an assumption which turns out to be consistent) one may approximate the factor /(pi, C)/(P2i C) m the 
integral (|B12j) as /(pi, C)/(P2, C) — / 2 (PiiC))> an d consider pi and p = p 2 — Pi as independent variables. Since in 
this regime f(p, £) is given by Eq. (jBlip the integral over pi and p 2 — P\ can be performed and the result takes is 



7T/ 2 



Do 



d( f d 2 q 



C 



q(z - C) 



exp 



iqp + 2k 2 



z-c 



dr) (g(m/k) - g(0)) 



(B15) 



Performing the angular part of the integral over q, expressing the pre-exponential factor as a derivative of the exponent, 
and changing the integration variable from ( to Z - ( wc finally obtain formula (I40|) : 



C(p) 



J 2 



ADgk 2 



C-z 



/d 
qdqJ {qp)—exp 



drj(l-g 



- (1 



(B16) 



where t 1 = fc^^O), while g(rf) = g±(v)/g±(0)- 
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